/*


   12/05/2010   Jacobi obtenido de 

	http://www.openmp.org/samples/jacobi.f


    Traducido por Kiko y Ruyman


*/


/*  program to solve a finite difference  */
/*  discretization of Helmholtz equation :   */
/*  (d2/dx2)u + (d2/dy2)u - alpha u = f  */
/*  using Jacobi iterative method.  */
/*   */
/*  Directives are used in this code to achieve paralleism.  */
/*  All do loops are parallized with 'static even' scheduling to  */
/*  maximize performance.  */
/*   */
/*  Input :  n - grid dimension in x direction  */
/*           m - grid dimension in y direction */
/*           alpha - Helmholtz constant (always greater than 0.0) */
/*           tol   - error tolerance for iterative solver */
/*           relax - Successice over relaxation parameter */
/*           mits  - Maximum iterations for iterative solver */
/*  */
/*  On output  */
/*        : u(n,m) - Dependent variable (solutions) */
/*        : f(n,m) - Right hand side function  */
/* ************************************************************ */

#ifndef LLC_TRANSLATION
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#include <math.h>
#endif
/* Constants */
#define N 4096
#define M 4096
#define ALPHA 0.00000005L
#define RELAX 0.01L
#define TOL 0.00001L
#define MITS 200



typedef struct timespec CLOCK_TYPE;

void chrono_gettime(CLOCK_TYPE *tv) {
        clock_gettime(CLOCK_MONOTONIC, tv);
}

double diftime(CLOCK_TYPE ch, CLOCK_TYPE ch2) {
   CLOCK_TYPE temp;
   if ((ch2.tv_nsec-ch.tv_nsec)<0) {
       temp.tv_sec = ch2.tv_sec-ch.tv_sec-1;
       temp.tv_nsec = 1000000000+ch2.tv_nsec-ch.tv_nsec;
   } else {
       temp.tv_sec = ch2.tv_sec-ch.tv_sec;
       temp.tv_nsec = ch2.tv_nsec-ch.tv_nsec;
   }
   return (double) temp.tv_sec + (double) (temp.tv_nsec*1e-9f);
}
CLOCK_TYPE start_total, end_total, start_kernel, end_kernel;


/* ***************************************************** */
/*  Initializes data  */
/*  Assumes exact solution is u(x,y) = (1-x^2)*(1-y^2) */
/*  */
/* ***************************************************** */
void initialize(int n, int m, double alpha, double *_dx, double *_dy,
	   double * u, double * f)
{
    int i, j, xx, yy;
    double dx = 2.0 / (n - 1);
    double dy = 2.0 / (m - 1);
/*  Initilize initial condition and RHS */

    //printf("*** Initialize  n: %d m: %d alpha : %g\n", n, m, alpha);
/* $omp parallel do  schedule(static)  */
/* $omp& shared(n,m,dx,dy,u,f,alpha) */
/* $omp& private(i,j,xx,yy) */
    for (i = 0; i < m; i++) {
	for (j = 0; j < n; j++) {
	    xx = -1.0 + dx * (double) (i);	/* -1 < x < 1 */
	    yy = -1.0 + dy * (double) (j);	/* -1 < y < 1 */
	    u[i * m + j] = 0.0;
	    f[i * m + j] = -alpha * (1.0 - xx * xx) * (1.0 - yy * yy) - 2.0 * (1.0 - xx * xx) - 2.0 * (1.0 - yy * yy);
	}
    }
/*  $omp end parallel do */

    //printf("*** Initialize  n: %d m: %d alpha : %g\n", n, m, alpha);
    *_dx = dx;
    *_dy = dy;
}


/* ***************************************************************** */
/*  Subroutine HelmholtzJ */
/*  Solves poisson equation on rectangular grid assuming :  */
/*  (1) Uniform discretization in each direction, and  */
/*  (2) Dirichlect boundary conditions  */
/*   */
/*  Jacobi method is used in this routine  */
/*  */
/*  Input : n,m   Number of grid points in the X/Y directions  */
/*          dx,dy Grid spacing in the X/Y directions  */
/*          alpha Helmholtz eqn. coefficient  */
/*          omega Relaxation factor  */
/*          f(n,m) Right hand side function  */
/*          u(n,m) Dependent variable/Solution */
/*          tol    Tolerance for iterative solver  */
/*          maxit  Maximum number of iterations  */
/*  */
/*  Output : u(n,m) - Solution  */
/* **************************************************************** */

void jacobi(int n, int m, double *_dx, double *_dy, double alpha, double omega,
       double *u, double *f, double tol, double maxit)
{
    int i, j, k;
    double error, resid, ax, ay, b;
    double * uold = malloc(sizeof(double) * n * m);

    double dx = *_dx;
    double dy = *_dy;


    /*  Initialize coefficients */
    ax = 1.0 / (dx * dx);	/*  X-direction coef */
    ay = 1.0 / (dy * dy);	/*   Y-direction coef */
      b = -2.0 / (dx * dx) - 2.0 / (dy * dy) - alpha;	/*  Central coeff   */

    error = 10.0 * tol;
    k = 1;
    chrono_gettime (&start_kernel);
#pragma acc kernels name("compute") copy(u[n*m], uold[n*m], error) copyin(f[n*m], omega, tol, n, m, ax, ay, b, alpha,i,j,resid)
     while (k < maxit && error > tol) { 
	error = 0.0;
{
    #pragma acc loop /* shared(omega,error,tol,n,m,ax,ay,b,alpha,uold,u,f)*/ private(i, j, resid) collapse(2)
	for (i = 0; i < m; i++)
	    for (j = 0; j < n; j++)
		uold[i * m + j] = u[i * m + j];
/*  Copy new solution into old */

    #pragma acc loop /*  shared(omega,error,tol,n,m,ax,ay,b,alpha,uold,u,f)*/ private(i, j, resid)  reduction(+ : error)
	for (i = 1; i < m - 1; i++) {
	    for (j = 1; j < n - 1; j++) {
		   resid = (ax * (uold[(i - 1) * m + j] + uold[(i + 1) * m + j])  /*      Evaluate residual  */
		    + ay * (uold[i * m + j - 1] + uold[(i * m) + j + 1])
		    + b * uold[i * m + j] - f[i * m + j]) / b;
		   u[i * m + j] = uold[i * m + j] - omega * resid;  /*  Update solution  */
		   error += resid * resid;                /*  Accumulate residual error */
	    }
    }

}
/* $omp end paralleldo  */
//   #pragma acc update device(error)
	/*  Error check  */
	k++;
        printf ("El error ahora vale %4.14lf\n", error);
	error = sqrt(error) / (double) (n * m);
    } /*  End iteration loop  */
    
    chrono_gettime (&end_kernel);
    //printf("Total Number of Iterations %d \n", k);
    //printf("Residual %g \n", error);
    *_dx = dx;
    *_dy = dy;

}



/* *********************************************************** */
/*  Checks error between numerical and exact solution          */
/* *********************************************************** */
void error_check(int n, int m, double alpha, double *_dx, double *_dy,
	    double * u, double * f)
{
    int i, j;
    double xx, yy, temp, error;

    double dx = 2.0 / (n - 1);
    double dy = 2.0 / (m - 1);
    error = 0.0;

/* $omp parallel do schedule(static) */
/* $omp& shared(n,m,dx,dy,u,error) */
/* $omp& private(i,j,xx,yy,temp) */
/* $omp& reduction(+:error) */
    for (i = 0; i < m; i++)
	for (j = 0; j < n; j++) {
	    xx = -1.0L + dx * (double) (i - 1);
	    yy = -1.0L + dy * (double) (j - 1);
	    temp = u[i * m + j] - (1.0L - xx * xx) * (1.0L - yy * yy);
	    error += temp * temp;
	}
    error = sqrt(error) / (double) (n * m);
    //printf("Solution Error : %g \n", error);
    *_dx = dx;
    *_dy = dy;
}

/* ************************************************************ */
/*  Subroutine driver ()  */
/*  This is where the arrays are allocated and initialzed.  */
/*  */
/*  Working varaibles/arrays  */
/*      dx  - grid spacing in x direction  */
/*      dy  - grid spacing in y direction  */
/* ************************************************************ */
void driver(int m, int n, int iters)
{
    double dx, dy;
   double * u;
   double * f;

   u = malloc(sizeof(double) * n * m);
   f = malloc(sizeof(double) * n * m);

    initialize(n, m, ALPHA, &dx, &dy, u, f);
    jacobi(n, m,  &dx, &dy, ALPHA, RELAX, u, f, TOL, iters);
    error_check(n, m, ALPHA, &dx, &dy, u, f);

}



int main(int argc, char *argv[])
{
   if (argc != 4) {
      printf("*** Incorrect number of parameters, expected progname <n> <m> <iter> \n");
      return EXIT_FAILURE;
   }
   int n = atoi(argv[1]);
   int m = atoi(argv[2]);
   int iters = atoi(argv[3]);
   chrono_gettime (&start_total); 
    driver(n, m, iters);
   chrono_gettime(&end_total);
    double kernel = diftime(start_kernel, end_kernel);
    double total = diftime(start_total, end_total);
    printf ("Tiempo kernel: %lf s\n", kernel);
    printf ("Tiempo total: %lf s\n", total);


   return 0;
}


